Quantum spin Hall effect in three dimensional materials: 
Lattice computation of Z2 topological invariants and its application to Bi and Sb 
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We derive an efficient formula for Z2 topological invariants characterizing the quantum spin Hall 
effect. It is defined in a lattice Brillouin zone, which enables us to implement numerical calculations 
for realistic models even in three dimensions. Based on this, we study the quantum spin Hall effect 
in Bi and Sb in quasi-two and three dimensions using a tight-binding model. 
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Quantum spin Hall (QSH) effect [J, i, S i, i] has 
been attracting much current interest as a new device 
of spintronics H, 0, S, Hi. It is a topological insulator 



[lOl I 111 . Il2l | analogous to the quantum Hall (QH) effect, 
but it is realized in time-reversal (T) invariant systems. 
While QH states are specified by Chern numbers 3, 14 1 , 
^SH states are characterized by Z2 topological numbers 

Graphene has been expected to be in the QSH phase 
P, However, recent calculations have suggested that 
the spin-orbit coupling in graphene is too small to reveal 
the QSH effect experimentally 3, IB- Recently, it has 
been pointed out that Bi thin film is another plausible 
material for QSH effect ■ Also by the idea of adiabatic 
deformation of the diamond lattice, it has been conjec- 
tured that Bi in three dimensions (3D) is in a topological 
phase [18]. 

While systems in two dimensions (2D) are character- 
ized by a single Z2 topological invariant, four independent 
Z2 invariants are needed in 3D [l^, [l^ 0] • This makes 
it difficult to investigate realistic models, in which com- 
plicated many-band structure is involved. Therefore, for 
the direct study of Bi in 3D as well as for the search for 
other materials, to establish a simple and efficient com- 
putational method of Z2 invariants in 3D is an urgent 
issue to be resolved. 

In this paper, we present a method of computing Z2 
invariants based on the formula derived by Fu and Kane 
2ll | together with the recent development of computing 



Chern numbers in a lattice Brillouin zone [23, l23|, |24|. 
This method is simple enough to compute Z2 invariants 
even for realistic 3D systems. Based on this, we study a 
tight-binding model for Bi and Sb. 

First, we derive a lattice version of the Fu-Kane for- 



mula 2l|. To this end, we restrict our discussions, for 



simplicity, to systems in 2D, where a single Z2 invari- 
ant is relevant. Let T be the time- reversal transfor- 
mation T = ia^K, and assume that the Hamiltonian 
in the momentum space 7i{k) transforms under T as 
TH{k)T-^ = n{-k). Let V(fc) = (|l(fc)), • ■ ■ , |2A/(fc))) 



denote the 2M dimensional ground state multiplet of the 
Hamihonian: 7^(fc)|r^(fc)) = Er,{k)\n{k)) [ul, As- 
suming that the many-body energy gap is finite, we focus 
on topological invariants under the U(2M) transforma- 
tion 



il){k) ->i:{k)U{k), C/(fc)eU(2M). 



(1) 



As discussed the pfafRan defined by pik) = 

Yii'^\T'^) characterizes the topological phases of T in- 
variant systems. To be precise, the systems belong to 
topological insulator if the number of zeros of the pfafRan 
in half the Brillouin zone is 1 (mod 2), and belong to sim- 
ple insulator otherwise. This number has been referred 
to as Z2 invariant. It should be noted that under Eq. 
([T]), the pfafRan p{k) transforms as p{k) e~^'^'^^^p{k), 
where 4'{k) is the U(l) part of U(2Af) deRned through 
the relation e'-^C') = deiU{kV 

Recently, Fu and Kane [2l| have shown that the Z2 
invariant is expressed alternatively by 
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where B" = [— tt, tt] ® [— tt, 0] (See Fig. [T]), and where 
A and F is, respectively, the Berry gauge potential and 
associated field strength defined hy A = Tnjj'^dip and 
F — dA [Tl|, [lij . Notice that the gauge transformation 
^ yields A ^ A-\- idcj). This implies that the gauge of 
the Berry gauge potential can be fixed by the condition 
that the pfaffian p{k) is real positive. Note also that 
A{—k) = A{k) holds in this gauge. This is a kind of T 
constraint, as stressed by Fu and Kane 2lj. The zeros of 
p{k) thus serve as an obstruction of the gauge fixing [2^. 

In systems with breaking T symmetry like QH effect, 
such an obstruction gives in general a nontrivial Chern 
number. Contrary to this, in systems under considera- 
tion, the Chern number always vanishes due to T invari- 
ance. Even in this case, an obstruction occurs in the Bril- 
louin zone as long as the zeros of the pfaffian exist. Since 
the zeros occur at the time reversed pairs of points ±fc*. 
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FIG. 1: Left: The Brillouin zone, where shaded region de- 
notes B" . The thick lines indicate the boundary of B~ . The 
integration over Cj gives a winding number. Right: A lattice 
on the Brillouin zone. The sites in B~, Bt , and are de- 
noted by blue, red, and black circles, respectively. The shaded 
region denotes the plaquettes in B~ . 



the vortices around these pairs are opposite and there- 
fore cancel each other, giving vanishing Chern number. 
Nevertheless the sum of vorticities in half the Brillouin 
zone, e.g., in B~ , just gives the number of zeros of p{k) 
(mod 2). Imagine, for example, that two zeros exist in 
B~ . Since they are generically first order zeros, the wind- 
ing numbers are ±1. Then their sum is restricted to ±2 
and 0, which can be denoted as "0 mod 2" . It thus turns 
out that the Fu-Kane formula ^ counts the vorticities 
in half the Brillouin zone. So far we have discussed in 
a specific gauge, but in any other gauge, D changes by 
2, provided that the gauge keeps T constraint. There- 
fore, D is indeed a Z2 topological invariant. It has also a 
topological stability against small perturbation as long as 
the many-body gap is finite. As we will show below, this 
expression for Z2 invariant is convenient for numerical 
computations. 

Define a lattice on the Brillouin zone, 



ke = 7r(ji/7Vi, j2/iV2), = -N^, ■■■ ,N^. (3) 

The sites labeled by k( are divided into three sets, 
and T invariant sites B^ denoted by red, blue and black 
circles in Fig. [1] respectively. Here, T invariant sites 
are specified by the property that TH{ke)T^^ — H{ke). 
As a T constraint, we choose the states at —kg e B^ as 
their Kramers doublets at kg E B^ . Suppose that at kg 
the spectrum is arranged as en{kf ) < En+iikg). Then the 
states at —kg can be constrained as 



\n{^kg)) ^T\n{kg)), for kg € B; 



(4) 



On the other hand, both of the Kramers doublets are in- 
cluded in : The spectrum in this set can be arranged in 
general as e2n-i{k) = £2n{k) < e2n+i(fc) • • • • Therefore, 
we enforce the constraint 

\2n{kg)) =T\2n-l{kg)), for kg e B^ (5) 

With these constrained states, we define a link variable 
=AA^"'(fc€)deti^t(fc^)^(^^ + ^)^ (6) 



where Af^ ^{ke) = \ det 1p^kg)■^p{ki + fl)\, and associated 
field strength through a plaquette variable 

Fu{kg) - lnUiikg)U2ikg + i)U^\kg + 2)U^\kg), (7) 

where F12 is defined within the branch Fi2/i G (— 7r,7r). 

The sum of i^i2 over B~ can be written as a similar 
formula to Eq. ([2|). To see this, it is convenient to de- 
fine a gauge potential via Af^{kg) = hiUfj,{kg) also in the 
branch A^{kg)/i € (— 7r,7r). Then the field strength can 
be rewritten as 

Fi2{kg) = AiA2{kg) - A2Ai{kg) + 2TTini2{kg), (8) 

where integral field ni2{kg) has been introduced so as to 
match the branches of both sides 27, 2§|. Thus, we 
reach 
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ni2{kg), (9) 



where the sums of F12 and of ni2 are over the plaquettes 
in the shaded region denoted by B~ in Fig. [1] The sum 
of Afj^ is over the links of the boundary of B^ specified 
by thick lines in Fig. [T] Therefore, a lattice version of D 
is 
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This formula for the Z2 invariants is one of the main re- 
sults of this paper. Indeed this formula has the following 
desired properties. Firstly, it is strictly integral. Sec- 
ondly, though the ground state multiplet can be mixed 
by Eq. ([1]), it is SU(2M) invariant. Finally, it changes by 
2 under the remaining U(l) transformation, and hence, 
it is Z2 invariant. The last property will be proved else- 
where, though it is not difficult. 

In 3D, it has been shown that the phases of T invariant 
systems are classified by four independent Z2 invariants 
flil . [TqI . [2oI |. To compute them, let us define six two- 
dimensional tori, according to Moore and Balents [3]. 
For example, fix the third momentum to /C3 = or tt, 
then we have two tori spanned by ki and ^2 which we 
denote Zq and torus, respectively. Applying the pre- 
vious techniques, we can compute two Z2 invariants Dl 
which are referred to as and z^^. In the same way, we 
have six invariants xq, x^, uq, y-^, zq, and living on 
six tori Xq, X^^, Yq, F^, Zq, and Z^r, respectively. There 
are two constraints, however: xqXtt — y^^y-K — -Zo-^tt (mod 
2), and therefore, four invariants among six are indepen- 
dent [3]. According to Fu et al. [1^, we choose them 
as i/Q = xqx^, vi — Xtt, V2 — y-K, and z/3 = z^r, and 
denote them as vq\ {viV2V2,). As is known in the QH ef- 
fects, non-trivial structures of topological ordered states 
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are hidden in the bulk and play physical roles near the 
boundaries as chracteristic edge states [2^. Based on the 
principle, by investigating the relationship between the 
Z2 invariants and surface states, Fu et al. have clarified 
that there are basically three phases; simple band insula- 
tor, weak topological insulator (WTI) which is topolog- 
ical but weak against disorder, and more robust strong 
topological insulator (STI) jlsi]. 



FIG. 2: The n-field configurations for Bi (left two panels) and 
Sb (right two panels) computed by the gauge that the pfaffian 
is real positive. The first and third (second and fourth) show 
the configurations on the Zq (Z-n) torus. The white and black 
circles denote ni2 = 1 and —1, respectively, while the blank 
denotes 0. These read zo = 2 and = for Bi, and zo = 
and z-n = 1 for Sb. 

Recently, Murakami T?| has pointed out the possibil- 
ity of QSH effect in Bi. Though Bi is a semimetal, the 
valence band and conduction band keep the direct gap 
throughout the Brillouin zone. Fu et al. have studied 
solvable tight-binding models with the diamond struc- 
ture, and predicted that the valence band of Bi is char- 
acterized by the WTI phase specified 0;(111), based on 
the observation that the structure of Bi can be viewed 
as an adiabatically distorted cubic lattice toward the di- 
amond lattice. However, since a realistic tight-binding 
model including s and p orbitals with nearest neighbor, 
second neighbor, and third neighbor hoppings has indeed 
complicated band structure, we calculate the Z2 invari- 
ants directly for heavy group V elements. 
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FIG. 3: Phase diagrams of Bi (left) and Sb (right) as func- 
tions of Vp-pa and Vnp-n . Other parameters are the same 
as those in Ref. [29( ]. The colored regions denote differ- 
ent phases: red for 1;(111), yellow for I;(000), lightblue for 
0;(1I1), and blue for 0;(000). The white points indicate the 
position of the parameters for Bi and Sb in Ref. poj . 

In what follows, we investi gate real materials by the 
tight-binding models in Ref. [29[. We first discuss the 



phase of Bi which is attracting much interest. We show 
in Fig. [2] examples of the n-field configuration computed 
for Bi. Though these calculations are for rather coarse 
10 X 10 lattice, we have checked that finer ones indeed 
reproduce the same Z2 invariants and our formula is ac- 
tually convergent even in this mesh size. For Bi, these 
figures tell that zq = z^r = mod 2. The other Z2 in- 
variants are also mod 2, and it turns out that the va- 
lence bands of Bi in 3D is specified by 0;(000) phase. 
This result is contradictory to the conjecture by Fu et 
al. mentioned-above. This suggests that along adiabatic 
distortion of the lattice, some topological changes should 
occur. Indeed, a slight change of the Slater-Koster pa- 
rameters can lead_to different phases. Among adjustable 
14 parameters [29], crucial ones may be Vppa and Ypp^^, 
nearest neighbor hopping parameters between p orbitals. 
We show in Fig. [3] the phase diagram of Bi to discuss 
the stability of the phase. This diagram tells that Bi 
is located in a small island of 0;(000) phase embedded 
in 1;(111) phase. We also understand this feature from 
a small direct gap of Bi, 12 meV, at the L point [2^. 
With varying the parameters, this gap soon closes and 
the phase of Bi changes from trivial phase into STI phase. 
We conclude that Bi in 3D dose not show the QSH ef- 
fect, though it locates quite near the phase boundary 
with STI. 
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FIG. 4: Band structure of Bi (upper) and Sb (lower). Bi: 
The blue lines are energies at / = 0.5 and the red-dotted 
lines are at / = 0.1. Sb: The blue lines at / = 0.7 and the 
red-dotted lines at / = 0.1. There are other bands around 
— 10 eV associated mainly with s orbitals which are omitted 
in these figures. 
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On the other hand, decreasing the thickness, a 
semimetal-semiconductor transition occurs, and Mu- 
rakami has suggested that Bi thin film would be in QSH 
phase [17[. To study the quasi 2D systems, and also to 



clarify the discrepancy between our result and the con- 
jecture by Fu et ai, we next investigate the effects of 
dimensionality on the present model. We replace the 
Slater-Koster parameters (a = ssa, spir, ppa, ppir) for 
second neighbor hopping in Ref. [1^ by fV^ with a uni- 
form factor /. This factor / can effectively control the 
dimensionality into [111] direction, interpolating between 
Murakami's model at / = and the 3D model at / = 1. 

We show in Fig. |4] the band structure of Bi with 
/ = 0.5 and 0.1. At / = 1 (See Ref. f29|), the over- 
lap energy (indirect gap) is Ai? = — 12 meV. With de- 
creasing /, a semimetal-semiconductor transition occurs 
at / ~ 0.99. Fig. 2] confirms that Bi is indeed a semi- 
conductor at / = 0.5 and 0.1 whose overlap energy is 
AE = 56 and 90 meV, respectively. Near / = 0.993, the 




In Fig. [5l we show the phase diagram of Sb for / = 
0.1. However, it should be stressed that with appropriate 
thickness, 0.54 < / < 0.89, Sb is a semiconductor in STI 
phase and hence should show QSH effect. 

Finally, we comment on the relationship between the 
method presented in this paper and the previous one in 
Ref. 2j|. While in the present calculation link variables 
are defined with respect to the momentum, the previous 
calculation has been implemented with respect to twist 
angles by imposing a spin-dependent twisted boundary 
condition. For systems with appropriate strength of spin- 
orbit coupling, the present computation is more efhcient, 
but for systems with very small spin-orbit coupling as 
well as with inversion symmetry, the previous computa- 
tion by the use of the twisted boundary condition gives 
more reliable results. In this sense, both methods are 
complementary to each other. Details will be published 
elsewhere. 

We also mention that recently Fu and Kane 30] have 
reached the similar conclusion of the phases of Bi and Sb 
by making the use of inversion symmetry of the system. 
We stress here that our method can apply to any systems, 
even without inversion symmetry. 

This work was supported in part by Grant-in- Aid for 
Scientific Research (Grant No. 17540347, No. 18540365) 
from JSPS and on Priority Areas (Grant No. 18043007) 
from MEXT. YH was also supported in part by the Sum- 
itomo foundation. 



FIG. 5: Phase diagrams of Bi (upper panel) and Sb (lower 
panel) for / = 0.1 as functions of Vppa and Vpp,r. 

phase of Bi changes from 0;(000) into 1;(111). With fur- 
ther decreasing / and enhancing the two-dimensionality, 
topological change occurs again near / ~ 0.223, and the 
system becomes 0;(111), as shown in Fig. [5] This is just 
the phase predicted by Fu et ai [18]. Therefore, we sug- 
gest that the adiabatic distortion of the diamond lattice 
leads to Bi thin film, and along the change of the dimen- 
sionality, the adiabatic distortion does not work, giving 
rise to gap-closing and resultant topological changes. We 
also conjecture that STI phase is very stable along the 
change of /, and could be observed by experiments. 

Sb is also a semimetal with a larger gap at the L point 
[2^. We show in Fig. [3]the phase diagram of Sb. It turns 
out that Sb belongs to the 1;(111) phase even in 3D. Its 
location is far from the phase boundary with 0;(000) and 
therefore, it is rather stable. We show in Fig. 2] the 
band structure of Sb. With decreasing /, a semimetal- 
semiconductor transition also occurs at / ^ 0.89. Along 
the change of /, topological change occurs once: Near 
/ ^ 0.54, the phase changes from 1;(111) into 0;(000), 
and no 0;(111) phase is observed throughout. Therefore, 
the phase of Sb thin film is 0;(000), different from Bi. 
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